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ABSTRACT: Lattice quantum chromodynamics provides first principles calculations for 
hadrons containing heavy quarks — charm and bottom quarks. Their mass spectra, decay 
rates, and some hadronic matrix elements can be calculated on the lattice in a model indepen- 
dent manner. In this review, we introduce the effective theories that treat heavy quarks on the 
lattice. We summarize results on the heavy quarkonium spectrum, which verify the validity of 
the effective theory approach. We then discuss applications to B physics, which is the main 
target of the lattice theory of heavy quarks. We review progress in lattice calculations of the B 
meson decay constant, the B parameter, semi-leptonic decay form factors, and other important 
quantities. 
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1 Introduction 

Quantum chromodynamics (QCD) is the fundamental theory of the strong inter- 
action. It is the SU(3) gauge field theory, which describes the interaction among 
quarks carrying three (R, G, and B) internal degrees of freedom, conveniently 
called "color." Because of the non-Abelian nature of gauge symmetry, the gluon 
(the gauge field in QCD) interacts with itself, and as a result the strong coupling 
constant as becomes large as the momentum of the exchanged gluon decreases. 
Therefore, the calculation of low energy (approximately several hundred MeV) 
properties of hadrons, such as their masses and interactions, is a challenging 
problem that requires nonperturbative methods. 

A regularized formulation of QCD on a four-dimensional hypercubic lattice, 
called lattice QCD (1), enables us to calculate such hadron properties in the 
strong coupling regime. Because the interaction is highly nonlinear and the num- 
ber of degrees of freedom is huge (proportional to the space-time volume), the 
lattice calculation is computationally so demanding that several dedicated su- 
percomputers have been developed around the world. The application of lattice 
QCD covers a wide area including the light hadron mass spectrum, finite temper- 
ature phase transition, and weak-interaction matrix elements — as well as heavy 
quark physics, which is the subject of this article. 

Among the six flavors of quarks in the Standard Model, charm (c), bottom (6), 
and top (t) quarks are heavy compared to the hadronic energy scale. The top 
quark decays very rapidly to a real W boson and bottom quark, and its QCD 
interaction can be treated perturbatively. Therefore, we consider the physics of 
charm and bottom quarks. 

The main motivation for charm and bottom quark physics is to study the 
flavor structure of the Standard Model including the mechanism that induces CP 
violation. Detailed study of the decays of charmed and bottom mesons provides 
pieces of information to precisely determine the Cabibbo-Kobayashi-Maskawa 
(CKM) matrix elements, which are fundamental parameters in the Standard 
Model. Furthermore, the loop effect of physics beyond the Standard Model may 
be probed through flavor-changing neutral current (FCNC) processes, such as 
mixing or rare B decays, which are suppressed in the Standard 
Model by the Glashow-Iliopuolos-Maiani (GIM) mechanism (2). In such studies 
model-independent calculations of the hadronic matrix elements are crucial in 
order to precisely predict the Standard Model contributions. 

One example is the rate (or frequency) of the neutral B meson mixing, which 
is induced by AB = 2 box diagrams. Because of the GIM mechanism, the 
dominant contribution comes from the diagram involving the top quark, which is 
proportional to the CKM matrix element \Vtd\ squared. At the B meson energy 
scale, this interaction is effectively described by a local operator (9^^=2 = ^7^(1 — 
75)^67^(1 — 75)9, and its rate is proportional to a matrix element {B\0^^^^\B), 
which is parametrized by the B meson decay constant Jb and the S-parameter 
Bb- This matrix element represents the probability of flnding the heavy quark 
and the light anti-quark at the same spatial point in the B meson to annihilate 
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Figure 1: Typical momentum scales in the heavy- light (left) and heavy- heavy 
(right) mesons. 

it and to create a B meson from that point. This is a problem of bound state 
formation in QCD, which is non-perturbative. The lattice calculation offers the 
best tool to solve such problems, but one needs a special treatment of heavy 
quark fields on the lattice, since the Compton wavelength of the heavy quark, 
~ l/niQ, could be smaller than the typical lattice spacing a. 

The typical energy or momentum scale that governs the dynamics of quarks and 
gluons inside the usual light hadrons is the QCD scale Aqcdi which characterizes 
the energy scale where the QCD coupling Og becomes large. Because the heavy 
quark introduces a new energy scale to the system, the QCD dynamics inside 
heavy hadrons is quite different from that of light hadrons. 

In the heavy-light meson, which is a bound state composed of a heavy quark 
Q and light degrees of freedom (a light anti-quark q and gluons g), the motion 
of the heavy quark of mass tjiq is hardly affected by the light degrees of freedom 
with a typical momentum Aqcd, if ?71q » Aqcd- Thus, the heavy quark stays 
almost at rest at the center of the bound state, surrounded by the light degrees of 
freedom, as shown in the left panel of Figure 1. The motion of the heavy quark is 
suppressed by Aqcd/"1q. This system is described by the Heavy Quark Effective 
Theory (HQET) (3-5), which is discussed in Section 2. 

The dynamics of quarkonium, which is a bound state of Q and Q, is governed 
by different energy scales. In the classical picture, the nonrelativistic kinetic en- 
ergy (p^)/2mQ and the potential energy — |as(i) have to be balanced, and the 
heavy quarks move around each other, in contrast to the heavy-light dynamics, as 
depicted in Figure 1 (right panel), (p) and (j,) denote the typical size of the rela- 
tive spatial momentum and distance between the two heavy quarks, respectively. 
The uncertainty relation (p) ~ implies {p) ~ agmq, which means that the 

typical velocity f ~ (p) / mq is of the order of the strong coupling constant as ■ 
Thus, there are three momentum scales in the heavy quarkonia: heavy quark 
mass rriQ, typical spatial momentum {p) ~ rnqv, and typical binding energy 
{p'^)/mQ ~ rnqv^. It is possible to construcnt an effective theory that explicitly 
separates these energy scales. This is known as non-relativistic QCD (NRQCD) 
(6-8). 

A number of lattice simulations for heavy quarks have been performed with 
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or without these effective theories. For the heavy-hght mesons, the simplest and 
best known appUcation is the calculation of the B meson decay constant. Other 
more complicated quantities include the i?-parameters and the semi-leptonic de- 
cay form factors. These quantities are not known from experimental data, and 
are necessary inputs to determine the CKM matrix elements. For the heavy- 
heavy mesons, on the other hand, the mass spectrum of cc and bb quarkonia is 
experimentally very well known, and thus provides a benchmark by which to test 
the validity of the lattice calculations. 

Until recently, because of the lack of computer power, most lattice simulations 
have been done in the quenched approximation, in which the effect of quark- 
antiquark pair creation (and annihilation) in the vacuum is neglected. This ap- 
proach could introduce uncontrolled systematic uncertainty, and hence the results 
should not be taken as a model independent prediction. However, the full QCD 
simulation without such an approximation has become realistic recently. We 
review some results from such calculations. 

This article is organized as follows. In Section 2 we introduce the effective 
theories used for heavy- light and heavy- heavy systems, and discuss their lattice 
formulations. Matching of the effective theory onto the continuum QCD is one 
of the key issues. Some results for quarkonia mass spectrum and decays are 
discussed in Section 3. Model independent calculation of the matrix elements 
required in B physics phenomenology is the main motivation for the study of 
heavy quark on the lattice. We review the current status in Section 4. Finally, 
the prospects for future lattice calculations are discussed in Section 5. 

We do not discuss the determination of the bottom quark mass using the 
lattice results for quarkoinum because it was well covered in a recent review by 
El-Khadra & Luke (9). Reviews presented at recent lattice conferences (10-13) 
and in other review articles (13, 14) are also useful sources of information. For 
textbooks on lattice gauge theory, see for example (16, 17). 

2 HQET and NRQCD 

The lattice spacing a in today's typical lattice simulations is about 0.1-0.2 fm. 
The available physical volume is then y ~ (2 fm)^ for the typical lattice size 
L/a = 16-24, which is numerically feasible. Because the Compton wavelength of 
charm and bottom quarks {1/mc — 0.13 fm and l/irih ~ 0.04 fm, respectively) 
are similar to or even smaller than the lattice spacing, discretization errors would 
seem to be out of control. Naively, the discretization error grows as (amg)", 
and the power n is 2 for the 0(a)-improved lattice fermion actions, such as the 
0(a)-improved Wilson fermion (18). 

However, the non-perturbative dynamics of QCD becomes important only in 
the low energy (~ Aqcd) regime. Therefore, if one can formulate the theory 
such that the energy scale of order mq is explicitly separated from the low en- 
ergy degrees of freedom, it would suffice to treat only the low energy part in the 
lattice simulation; the high energy part can be reliably treated in perturbation 
theory. Theoretically, such separation of different energy scales can be formu- 
lated in terms of the Operator Product Expansion (OPE) (19). Its well-known 
application is in deep inelastic scattering, for which the large energy scale comes 
the energy of collisions. Lagrangian-based formulations are more transparent in 
many applications, and for the heavy-light and heavy-heavy hadrons dedicated 
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formulations have been developed: Heavy Quark Effective Theory (HQET) for 
the heavy-light and nonrelativistic QCD (NRQCD) for the heavy-heavy hadrons. 



2.1 Continuum HQET and NRQCD 



Let us consider a heavy hadron in its rest frame. The momentum of the heavy 
quark inside the heavy hadron can be written as Pq = {mq + k^, k) by separating 
the heavy quark mass mq from other, smaller momentum. From the heavy 
quark field q one may extract the trivial dependence on the heavy quark mass 

as q = e"*"**?* ^ ^ ^ ' where the four-component spinor is divided into two two- 
component spinors Q and X. Then, the Dirac equation (ij^Df^ — mq)q = with 
the Dirac representation of the 7 matrices can be separated into two parts: 



iDoQ 
{2mq + iDo) X 



ia-DX, 
ia ■ D Q, 



(1) 
(2) 



where are the Pauli matrices. The second equation shows that the lower 
component X is smaller than the upper component Q by a factor 2mq. If we 
can neglect the small time-dependence of X, iDqX in Equation 2, and substitute 
X = ia ■ D/2mqQ into Equation 1, we arrive at the non-relativistic Shrodinger 
equation plus the Pauli term, 



iDoQ 



2nic 



+ 



gB 



'^rnq 



Q, 



(3) 



where g is the QCD coupling constant and -B* = ^e^^^F^^ is the magnetic part 
of the field-strength tensor of QCD. In Lagrangian form we may write this as 



eavy 



Q 



iDo + 



'^■rnq 



gB gE - gE D 



2m, 



Q 



Smi 



a ■ {iD X gE - gE X iD) {D 



2\2 



8mi 



8m'} 



+ 



Q, 



(4) 



recovering some of the higher order contributions. The ellipses represent ne- 
glected higher order terms. 

In the heavy-light hadron, the light degrees of freedom (light quarks and gluons) 
have momenta of order Aqcd- Exchange of spatial momenta with the heavy 
quark occurs through the 1/mq and higher order terms in Equation 4. Thus, the 
motion of the heavy quark is suppressed by powers of ^qcd /mq. In HQET (3- 
5) one treats the Lagrangian equation (Equation 4) as a systematic expansion in 
1/mq. At leading order, only the first term QUDqQ remains, which represents a 
non-moving heavy quark acting only as a static color source. This is often called 
the static approximation. The higher order terms can be included by operator 
insertions, that is, by expanding exp(i f d^xCheavy) in terms of l/mq. In the 
lattice simulations, it is more convenient to include the higher order terms as a 
part of the heavy quark propagator. 

In the heavy-heavy hadron, on the other hand, the momentum of the heavy 
(anti-)quark is determined by the balance between the potential energy and ki- 
netic energy, as already mentioned in Section 1. Specifically, in order to satisfy 
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Equation 3, {iDqQ) ~ {D"^ /IrnqQ) . Because the heavy quark potential is well 
described by the Coulomb form —Cpas/r, this implies 

(^) ~ {a.p) ~ ^. (5) 
r mq 

Then, the typical momentum is of the order mqas, or the typical velocity v = 
(p) /mq is of order . The separation of the energy scales is provided by the 
velocity v: (heavy quark mass mq) ^ (momentum mqv) ^ (binding energy 
mqv'^). 

In NRQCD (6-8) the counting of terms in the Lagrangian equation (Equa- 
tion 4) is done by the power of v. Each operator has a different power counting: 
Q ~ (mgi))"^/^, D ~ mqv, Dq ~ mqv^ , gE ~ vnqV^ , and gB ~ mqv^. As a 
result, the leading terms in Equation 4 are QHDqQ and {—D'^ /2mq) Q; the 
other terms are suppressed by a relative power v'^. Higher order terms can also 
be included in a systematic manner. 

The QCD coupling constant as depends on the energy scale, according to the 
renormalization group running. It should be evaluated at the momentum of the 
exchanged gluons, which in this case is the momentum scale of mqv ~ mqUs- By 
solving V ~ as{mqv) self-consistently, the typical velocity squared is estimated 
to be ~ 0.3 for charmonium (cc) and 0.1 for bottomonium (hb). Expansion in 
is therefore very effective for hb but marginal for cc. 



2.2 Lattice HQET and NRQCD 

Once the effective theories are introduced, the length scale to be treated on 
the lattice is not 1/mq, but 1/Aqcd for heavy-light or l/mqv for heavy-heavy 
systems. These are long enough to be approximated by a lattice of a ~ 0.1-0.2 fm. 
A rough sketch of the relevant length scales is shown in Figure 1. 

Discretization of the Lagrangian equation (Equation 4) is straightforward. 
First, one must perform the Wick rotation = ix^ to move on to the Euclidean 
space-time. Then, the covariant derivative is replaced by a finite difference 

A(,+)g(x) ^ U^{x)Q{x + /i) - Q{x), (6) 

where x + jl denotes a lattice point next to x in the /x-direction. U^{x) is the link 
variable, which represents the gauge field on the lattice: the relation to the con- 
tinuum vector potential A^{x) is approximately given by U^{x) = e-iC£>[iagA^{x)\. 
The discretized derivative. Equation 6, is constructed such that it is covariant 
under the SU(3) gauge transformation. One can also define 

A(,-)Q(x) ^ Q{x) - Ulix - (l)Q{x - fL), (7) 

A^i^^ = (A(+) + A(-))/2, and the lattice Laplacian A^^) = ^. aJ+^aS"^ 
One explicit form of the lattice action is 

S = Y^ ^) ^) - KtQit - 1, x)], (8) 
t,x 

where t = X4. The operator Kf is a kernel of time-evolution: 
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Here subscripts represent the time slice at which Hamiltonian operators such as 
(1 — aHQ/2n) act, and an integer n is introduced to suppress instabihty that 
appears in the evolution equation owing to unphysical momentum modes (7, 8). 
The leading order Hamiltonian Hq is given by 

A(2) 

^0 = , (10) 

and the higher order terms are 

a gB A'^^^ ■ gE - gE ■ A'^^^ a ■ {A'-^^ x gE - gE x A'^^^) 



5H 



2mQ 8mQ Sm^ 

The field strength tensors B and E are made from link variables connected on a 
plaquette. Details of their construction are given, for instance, in Reference (8). 
One can also define other heavy quark actions, which are equivalent to Equation 8 
up to discretization effects. 

Because the propagation of the heavy quark is simply expressed by the time- 
evolution kernel, Equation 9, numerical calculation of the heavy quark propagator 
is much faster than calculation of the light quark propagator, which requires some 
iterative solver, such as the conjugate gradient method. 

As in the continuum case, the difference between HQET and NRQCD is in 
the underlying dynamics. We can use the same lattice action, but the size of the 
contribution from each term in the Lagrangian differs between heavy-light and 
heavy-heavy hadrons. In other words, at a given accuracy (say Aqcd/^-q in the 
heavy-light, or v'^ in the heavy- heavy) , the terms to be included in the calculation 
are different. 

To recover the original four-component Dirac spinor -0/^ from the two-component 
fields, one must apply the Foldy-Wouthuysen-Tani (FWT) transformation 



7 . A(=^) A(2) S • gB i747 • gE 
2mQ StjIq StHq 



? ) > (12) 



with S-' = diagjo"-', a^}. It is necessary when we construct operators that involve 
both heavy and light quark fields, such as the heavy-light axial-vector current 

In the language of the functional integral to define the quantum field theory, 
the FWT transformation is nothing but a change of variables, and the physics 
described by the theories before and after the transformation is the same up to 
neglected higher order terms in Equation 12. There is a class of possible such 
transformations, but the choice in Equation 12 is special because it removes the 
lower components from the theory. This suggests that there are other effective 
theories for which the FWT transformation is partially applied. Such a formula- 
tion is known as the Fermilab action (20, 21). It uses the usual four-component 
spinors on the lattice, and the higher spatial derivatives are introduced according 
to the order counting of HQET or NRQCD. This formulation is useful, because 
at the leading order (of 1/?tiq or v"^) the action can be taken to be the same as 
that of light quark actions, i.e. the Wilson fermion with or without the 0{a)- 
improvement. Therefore, this formalism covers entire mass regions (from light to 
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heavy) with a single fermion action, although the achievable accuracy varies de- 
pending on the quark mass. In the heavy quark regime the accuracy is estimated 
in terms of the HQET or NRQCD counting rule. The discretization effect can also 
be estimated by interpreting the lattice effective action with the Symanzik effec- 
tive theory (21). To correctly incorporate the higher order corrections one must 
introduce higher derivative terms to the action just as in the NRQCD action. 
An equivalent formulation, but without the FWT transformation, was proposed 
later in a different context (22). 



2.3 Renormalization and continuum limit of the effective theories 

Although the HQET and NRQCD capture the physical picture of the heavy 
quark inside the heavy hadrons, the resulting field theory is a non-renormalizable 
effective theory. This means that an infinite number of terms are needed in order 
to eliminate ultraviolet divergences that appear in quantum loop calculations 
of the theory. In fact the Lagrangian equation (Equation 4) contains infinitely 
many terms at higher orders, and to renormalize them one needs an infinite 
number of renormalization conditions, i.e. input parameters. In practice, one 
is only interested in finite orders of l/mg (or f^) and truncates the Lagrangian 
at that order. Then, the number of renormalization conditions is finite and the 
calculation is feasible. This is how these effective theories work: at a certain 
order of the systematic expansion the number of parameters is still finite, and 
thus the predictive power remains. 

On the lattice with a finite lattice spacing no divergence appears, but we still 
need input parameters. The coefficients in front of the terms in Equation 11 
are correct only at tree level; if we include the quantum effects, they become 
renormalized. These input parameters are provided by calculating some physical 
amplitudes in the full relativistic theory. The coefficients are then determined 
such that the effective theory gives the same physical amplitude at a given order 
of 1/niQ (or u^). This procedure is called matching, and it is usually done in 
perturbation theory (perturbative matching). 

At leading order the necessary terms are Q^Q, Q^DoQ and Q^D^Q. Their 
renormalization corresponds to the energy shift, wave function renormalization, 
and mass renormalization, respectively. The perturbative matching of these pa- 
rameters onto continuum full theory was calculated for several specific definitions 
of lattice NRQCD actions (23-26) and for the Fermilab action (27), although the 
complete calculation of higher order terms, especially the term a ■ B, is still to 
be done. 

Similarly, the operators to be measured on the lattice must be matched onto 
the continuum full theory. For example, the heavy-light axial-vector current 
= 'fphl5lfM'>Ph which is related to the B meson leptonic decay constant, is 
renormalized as 



A4, = Za 



2mQ 2mQ 



(13) 



for the continuum heavy-light axial-vector current defined, for instance, using 

the MS scheme. The matching parameters Za, and are calculated at 
the one-loop level of perturbation theory. In the static limit (mq 00), the 
perturbative calculation of Za was carried out for both the Wilson (28-30) and 
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the 0(a)-improved Wilson (31, 32) fermions for the Hght quark, and even a 
non-perturbative calculation has been done recently (33). One-loop calculations 
are also available for the NRQCD (26, 34-36) and for the Fermilab (37) lattice 
actions. The matching of the four-quark operator Ol = ''Ah 75(1 — 7^t)'0/V'/i75(l ~ 
7/i)V'Z) which is needed for the the B^-B^ mixing matrix elements, has also been 
calculated (31, 38-42). 

The continuum limit of the lattice HQET/NRQCD is not trivial. Beyond the 
static approximation, the lattice action contains higher dimensional operators 
that induce power divergent terms such as l/(amQ)" in the radiative corrections. 
Because the divergences are unphysical and must be subtracted (or renormalized) , 
the perturbative matching becomes less accurate as a — > 0, even though the 
physical picture of the Aqcd/'"^q (or v"^) expansion does not change. Therefore, 
the lattice calculation is practically restricted to relatively large lattice spacings, 
such that aruQ ^ 1 is maintained. This means that the continuum limit a ^ 
cannot be reached within the effective theories beyond the leading order (static 
limit), and one must check the stability of the results in a limited region of lattice 
spacings to keep amq > 1. For the Fermilab action (20), on the other hand, the 
continuum limit can smoothly be reached, because it is reduced to the Wilson-like 
fermions as amq — > 0. 

In the dimensional regularization adopted in the continuum perturbation the- 
ory, the problem of power divergences in the OPE appears as the renormalon 
ambiguity (43); that is, the lowest order Wilson coefficient contains an ambiguity 
of order Aqcd/'^^q due to the poor convergence behavior of the perturbative ex- 
pansion. The renormalon ambiguity cancels with the higher order matrix element 
of the same order (Aqcd/'^t-q), and the physical quantity to be computed is free 
from the ambiguity. In the theory with an explicit cutoff, such as lattice gauge 
theory, the problem becomes the difficulty of subtracting unphysical power di- 
vergences using the perturbative expansion, as discussed above. In practice, the 
perturbative expansion converges well as far as one keeps amq > 1, but higher 
order perturbative calculations are needed to obtain better control of systematic 
errors (10, 13). 

One can avoid the whole problem by performing the matching calculation non- 
perturbatively, as recently proposed by the ALPHA collaboration (44). The 
matching of lattice HQET onto the relativistic action is feasible on a lattice of 
physically small volume, L ~ 0.2 fm, because one can work with a small lattice 
spacing for which the relativistic action can be safely used. The matching of 
HQET on larger volumes is then performed by comparing the lattices of size 
L and 2L recursively. At present, the method is applied to the leading-order 
HQET (static limit), but if the matching could be done in a similar manner 
for higher dimensional operators, it would provide a breakthrough to go beyond 
perturbative matching. 

2.4 Relativistic lattice actions used for heavy quarks 

In the usual lattice fermion actions, the quark mass is included as a perturba- 
tion to the massless limit. Because the theory is renormalizable, the number of 
relevant operators is limited and the transition toward the continuum limit is 
straightforward. 

However, in practical simulations with present computer resources, the lattice 
spacing a is not always small enough to satisfy amq <C 1 for the charm quark. 
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and it is absolutely too large for the bottom quark. The discretization error 
starts from O^amq) for the unimproved Wilson fermion, or from ©[(amg)^] for 
the 0(a)-improved actions. 

One can see how the discretization errors behave by taking the energy-momentum 
dispersion relation as an example. The non-relativistic effective Hamiltonian of 
the "relativistic" lattice fermion can be written as 

Q, (14) 

where the rest mass mi and the kinetic mass m2 can be different if the Lorentz 
invariance is lost owing to the discretization error. At tree level they become 

mio = ln(l+moa), 
1 _ 2 1 

moa(2 + TTT-oa) 1 + moa' 

for the 0(a)-improved and -unimproved Wilson fermions, respectively. The term 
mo denotes the quark mass appearing in the relativistic lattice action. Then the 
problem appears as a violation of the relation mo = mi = m2 of order (amg)^. 

One way to avoid the problem of this large amq error is to reinterpret the lattice 
action using HQET as done in the Fermilab formulation (20). But, if one does not 
employ this reinterpretation, one must restrict oneself to the region where amq 
is small enough that the remaining errors of order (amg)^ and higher are under 
control. This is possible for the charm quark if one takes the lattice cutoff much 
higher than 2 GeV and eliminates the leading error (amg)^ by an extrapolation to 
the continuum limit, although it is numerically rather demanding. Such analysis 
has been done only recently for the charm quark mass (45) and for the Ds meson 
decay constant (46) in the quenched approximation. 

To reach the bottom quark mass one must rely on an extrapolation in the 
heavy quark mass. It can be done by using the heavy quark scaling law, e.g. for 
heavy enough quark masses the heavy-light meson decay constant fp behaves as 
l/\/Mp up to power corrections. This method requires a careful analysis because 
the discretization error of order (amg)^ grows as mg becomes larger, and this 
growth could even be amplified by the extrapolation, producing uncontrolled 
systematic errors. An extrapolation to the continuum limit a — > must be done 
before the heavy quark extrapolation to reduce such a problem. A combined fit 
with the result in the static limit may be used to stabilize the extrapolation. 

Because the large energy scale of order mg flows only in the temporal direction 
in the rest frame, the discretization error appears from the temporal derivative 
of the lattice action. It is therefore possible to avoid the large discretization error 
by taking a temporal lattice spacing at much smaller than 1/mg while keeping 
the spatial lattice spacing a^. This is called the anisotropic lattice action (47). 
Although the radiative correction could reintroduce large a^mgOs contributions, 
in wihch case the virtue of the anisotropy is lost (22, 48), it is also possible to 
construct lattice actions free from this problem (49). 

3 Quarkonia 

Since the discovery of J/V'j experimental studies of heavy quarkonia have pro- 
duced enormously rich and precise data for their mass spectra and decay rates. 



.eavy 



nil 



2mc 



+ ••• 



(15) 
(16) 



Heavy Quarks on the Lattice 
GeV 
10.5 + 



11 



10.0-- 



.5-- 



-•-&- 



MeV 
20-- 

0- 

-20- 

-40- 



hh 



-i-f— Xb2 

Xbi 

-XbO 



'Pi 



Figure 2: Bottomonium spectrum from Reference (52). The spin dependent 
splittings of the IP states are shown separately in the right panel. Dashed lines 
indicate the experimental value. The lattice data are for Nj = {filled circles) 
and Nf = 2 {open circles). 



Theoretically, the quark potential model description has been successful in ex- 
plaining their properties, indicating that heavy quarkonium is a non-relativistic 
bound state. Therefore, the study of quarkonia gives a firm testing ground for 
lattice QCD, and the results may also be used to calibrate the basic parameter of 
the calculation, the lattice spacing a. Furthermore, the lattice QCD will be useful 
for studying the properties of non-standard states, such as qqg hybrids and qqqq 
multi-quark states, which may depend on the details of non-perturbative QCD 
dynamics. 



3.1 Mass spectrum 

Heavy quarkonium is a system for which one of the most accurate calculations 
can be done on the lattice. Indeed, Davies et al. demonstrated that the mass 
spectrum can be precisely calculated using lattice NRQCD (50-53). They used 
the lattice NRQCD action including the entire 0{v'^) terms as given in Equa- 
tions 10 and 11. The leading discretization errors are eliminated by adding the 
O(a^) improvement terms. The bottomonium spectrum obtained on the lattice 
with 1/a ~ 2.4 GeV is shown in Figure 2 (52, 53). The masses of low- lying radial 
and orbitally excited states are precisely computed in both quenched {Nj = 0) 
and unquenched {Nf = 2) QCD. It is remarkable that the 1S-2S and IS-IP split- 
tings become more consistent with the experimental value if dynamical quarks 
are introduced. 

With the NRQCD formulation one can investigate the effect of each term in- 
cluded in the action by comparing the results with and without that contribution. 
Furthermore, because the bottomonium system is mostly sensitive to short dis- 
tance physics, our experience from quark models and perturbation theory can be 
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used to estimate the systematic uncertainty due to truncated higher order terms 
in and in a. 

In the non-relativistic expansion the leading terms that are not included in 
the action employed by Davies et al. are those of 0{v^). Because the typical 
v'^ is about 0.1 for the bottomonium system, the relative error to the leading 
terms (~ v'^) is naively O(u^) ~ 1%. An estimate using the quark potential 
model suggests even smaller error (54). The leading discretization errors are the 
radiative correction to the O(a^) improvement term and the O(a^) corrections. 
These are 0{as{TT/a)a'^p'^) and O(a^p^), respectively, for typical spatial momen- 
tum p = rnqv. For the lattice spacings between 0.05 fm and 0.15 fm the size of 
their corrections is estimated to be at most a few percent, and lattice data at 
three lattice spacings support such small discretization effects (53). 

For the spin-dependent splittings, the relative error is an order of magnitude 
larger because the leading contribution itself is O(u^). Indeed, Manke et al. 
studied the effect of the tree-level 0{v^) corrections in quenched NRQCD and 
found that the spin-dependent splittings are affected by 10%-20% when 0{v^) 
terms are added to the action, whereas the spin-averaged splittings are essentially 
unaltered, supporting the above naive expectation (55). 

As shown in Figure 2 the mass spectrum significantly deviates from experiment 
unless the effect of dynamical quarks is included. The sea quark effects have also 
been studied (56, 57) and found to be sizable. The recent calculation by the 
HPQCD and UKQCD collaborations on the gauge ensembles with 2+1 flavor 
sea quarks produced by the MILC collaboration indicates that both the spin- 
averaged splittings and the hyperfine splittings are in excellent agreement with 
experiment (58, 59). 

Because the bottomonium spectrum can be precisely calculated with good 
control of the systematic uncertainty, it provides a reliable input for the lattice 
scale 1/a of the given lattice. Then, other physical quantities can be calculated in 
the unit of the bottomonium mass splitting, such as the IS-IP splitting. Among 
other quantities, the relation between and the lattice scale 1/a can be obtained 
rather accurately. Using two-loop perturbation theory to relate on the lattice 
with that in the continuum renormalization scheme, e.g., the MS scheme, one 

(5) 

can determine ajjg{Mz) using the low energy hadronic scale as input. In a 
recent simulation with 2+1 flavors of dynamical quarks Davies et al. obtained 
ajj^{Mz) = 0.121(3) (58), which agrees with other experimental determinations 
(60). 

The charmonium spectrum is much harder to calculate precisely using the 
NRQCD formalism, because of larger relativistic corrections (w^ ~ 0.3). With 
the same NRQCD action used by Davies et al. the higher order contribution is 
expected to be about 10% for spin-averaged and 30% for spin-dependent split- 
tings. In fact the next-to-leading order 0{v^) terms produce a sizable effect on 
the spectrum (61), and the prescription to (approximately) incorporate radiative 
corrections strongly affects the hyperfine splitting (62). Therefore, relativistic lat- 
tice actions have been used for the charm quark in many studies on anisotropic 
(63, 64) and isotropic (65) lattices. Extrapolation toward the continuum limit is 
essential to control the ©[(amg)^] errors in such calculations, and the unquenched 
simulation remains to be intractable. 
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3.2 Matrix elements 

After successful reproduction of the quarkonium spectra, the next interesting 
problem would be to compute the quarkonium decay matrix elements in lattice 
QCD. Bodwin et al. (66) have shown that the quarkonium decay amplitude can 
be factorized into a product of long-distance matrix elements of four-fermion op- 
erators and short-distance coefficients that are calculable in perturbation theory. 
For instance, the S- and P-wave decay rates are written as 

rC'+^Sj^X) = giC'+'Sj)2lmhe^+'Sj)/ml 

+ T^e'+'SJ)2lmg^e^+^SJ)/mt, (17) 

rC^+^Pj^X) = Wi(2^+iPj)2Im/i(2-+iPj)/m^ 

+ nse'+'PJ)2lmg^e^+'PJ)/mt, (18) 

where / and g are the short-distance coefficients and the quarkonium states are 
indicated by '^^'^^Lj with spin s, orbital angular momentum L, and total angular 
momentum J. The matrix elements are defined in terms of the bottom quark 
field and anti-quark field x in the NRQCD formalism as 

Gi = eSo\^PhxH\'So), (19) 

^1 = eSo\i^^xxH-'^D)^i^\'So), (20) 

ni = i'PM^'-Dxx^'-Di^l'Pi), (21) 

Hs = (^5o|^tr'^xx^T>|i5o), (22) 

where the derivative D is defined as x^ Dil) = x^{DiIj) — (DxVi^- 

The lattice calculation of these matrix elements has been carried out using the 
O(f^) lattice NRQCD action for both quenched (67) and unquenched (68) QCD. 
In order to obtain the matrix element in the continuum NRQCD, the lattice 
operators are matched to those in the continuum at the one-loop order of pertur- 
bation theory. Then, the results of Qi can be compared to the phenomenological 
estimate obtained from the leptonic decay width 

r(T ^ e+e-) ^ (l-'-^)g,=5.0± 0.2 GeV. (23) 

3m^ \ 3tt J 

Bodwin et al. found that the matrix element Qi in quenched QCD is 40% smaller 
than the phenomenological estimate, whereas the unquenched result is consistent 
with the phenomenological value within the error of 10%. Although further 
improvements are required, it is encouraging that one finds good agreement at 
the present accuracy of O(f^) ~ 10 % level. 

3.3 Exotics 

Most of the known excited states of bb and cc mesons are classified in terms of 
the quark model; that is, their J^'^ is explained by the standard assignment 
P = (—1)-^''"^ and C = (— l)-^"*""^ where 5 = 0, 1 is the total spin of q and q and 
L = 0, 1, 2, ... is the orbital angular momentum. In QCD, however, there could 
be other excited states, which may be explained by a gluon excitation qqg or a 
four-quark state qqqq. An experimental candidate of the latter, called X(3872), 
was recently found in the charmonium system by the BELLE collaboration (69). 
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Figure 3: Heavy quark potential including gluon excitations. The plot is from 
(70, 71). 
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Figure 4: Charmonium spectrum from Reference (78). The rightmost three points 
correspond to the hybrid states. 



The quark-anti-quark potential has been calculated on the lattice, including 
the gluon excitations (70, 72, 73). Figure 3 shows such a potential with and 
without first gluonic excitation in the quenched approximation. Starting from the 
measured potential one can calculate the levels of hybrid mesons by employing 
some assumptions, such as the Born-Oppenheimer approximation (74). 

It is also possible to directly calculate the hybrid spectrum in lattice QCD by 
choosing interpolating operators with appropriate quantum numbers. To obtain 
a good statistical signal in the lattice calculation, the anisotropic lattice is useful, 
and such a lattice calculation has been done (74-78). The cc spectrum obtained 
by (78) is shown in Figure 4. 

4 B physics 

The strongest motivation for the lattice simulation of heavy quarks comes from B 
physics. One of the goals of the B factory experiments is to precisely determine 
the CKM matrix elements, which are fundamental parameters of the Standard 
Model. In particular, the two smallest elements V^b and Vtd are interesting be- 
cause they describe the mixing between first and third generations and are the 
source of CP violation. Furthermore, by studying a variety of flavor-changing- 
neutral-current (FCNC) processes, which are suppressed in the Standard Model 
by the GIM mechanism, one may get insight into the physics beyond the Standard 
Model entering into the loop diagrams. 

The bottom quark decay inside the B meson is always accompanied by a gluon 
and quark cloud, which makes the extraction of fundamental parameters from 
experimental data difficult. Therefore, a model independent calculation of QCD 
effects by a theoretical method, such as lattice QCD, is essential. 

Owing to the practical limitation of computer power, the momentum scale that 
lattice QCD can treat is limited to < 1 GeV, which is why the effective theories as 
discussed in Section 2 are developed. Computer power also limits the momentum 
of initial and final state hadrons to be simulated on the lattice to < 1 GeV, 
whereas in realistic B decays the momentum of the daughter particles can be 
as large as ~ mB/2. This means that only a limited range of kinematical space 
is covered by lattice QCD. Another important limitation of lattice calculation 
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is that one cannot treat the multi-hadron final state on the lattice, because the 
relation between the correlator in the Euclidean space and the amplitude in the 
Minkowski space-time is non-trivial when more than one particles are involved. 
Therefore, the direct calculation of the decay amplitude of B ^ tttt, for instance, 
is far beyond the reach of present lattice calculation. 

In the following section, we review the lattice calculation of several important 
B meson matrix elements. 

4.1 B^-B^ Mixing 

The mass difference (or oscillation frequency) between two neutral B mesons 
AMb^ [q = d oy s) is given by 

AMs, = ^r^BMBjlBB,ml,So{mi/rn^^)\Vt,Vt,\\ (24) 

where Sq{x) is the Inami-Lim function (79) originating from the box diagram and 
rjB is a short distance QCD correction calculated in perturbative QCD. The non- 
perturbative quantity to be computed on the lattice is the combination, f^^Bsg- 
The leptonic decay constant fs^ is defined through 

ifB.Pf. = {0\A^\Bg{p)), (25) 

with a heavy-light axial-vector current = ^757^^6. The B parameter Bsg is 
defined as 

(S0|0^«=2(;.)|i?0) ^ iBsMfsMk (26) 
with the AB = 2 operator 

0^5=2 ^ -^^(^ _ ^^)bq^^{l - 75)6. (27) 

This operator has an anomalous dimension, and as a result the B parameter is 
scale dependent, Bs^ilj)- For convenience a scale-invariant B parameter, Bb^ = 
5B,a^(/z)-6/23(]^ + has/Air), with J5 = 5165/3174, is often quoted. 

Experimentally, the mass difference AMb is already known very precisely for 
the B^ meson {AMb = 0.502 ± 0.007 ps"^ (80)), and the main problem in 
extracting the CKM element \Vtd\ is in the theoretical calculation of J^Bb- For 
the Bs meson, to date only a lower bound AMg > 14.4 ps"^ (95% CL) (80) 
has been obtained. The leptonic decay B~^ — > li^i is hard to measure at existing 
experimental facilities because of its small branching fraction, although there is 
a good chance that it will be measured and thus a combination will be 

determined at future higher luminosity B factories. 

Because the B^i and Bg mesons differ only in the valence light quark mass as 
far as QCD is concerned, one can expect that the theoretical uncertainty largely 
cancels in the ratio 

AMb, \Vt,\^^ ' ^ ' 

up to the SU(3) breaking effect. Once the Bg-Bg mixing rate is experimentally 
measured with good precision, it will help to determine the most interesting CKM 
matrix element | Vtd \ (assuming the CKM unitarity relation \ Vts\ = \ Vcb \ ) ■ 
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The width difference AF^ of Bg mesons is induced by the final states into 
which both Bg and Bg mesons can decay. (The width difference of in the 
Standard Model is too small to be observed.) The main contribution comes from 
6s — > cc — > bs transitions at the quark level. Because large momentum of order 
mb/2 flows into final state particles, the amplitude can be expressed in terms of 
the heavy quark expansion (81, 82) and the leading order is represented by matrix 
elements of local AB = 2 four-quark operators Ol = S7^(l — 75)657^(1 — 75)6 
and Os = s{l — 75)65(1 — 75)6. The corresponding matrix elements are Bb and 
Bs, defined such that they are unity in the vacuum saturation approximation. 

4-2 B^-B^ Mixing: Quenched Lattice Results 

In the quenched approximation, the B meson decay constant has been studied 
extensively. Recent calculations using the Fermilab action (83-85), the NRQCD 
action (26, 87-89), and the extrapolation method (90-92) show good agreement 
within an error of 0(15%). In the recent lattice conferences, world averages read 
Jb = 173 ± 23 MeV, fs, = 200 ± 20 MeV, and fsJfBa = 1-15 ± 0.03 (11, 12). 
The statistical error of the Monte Carlo simulation is a minor part of the error 
given above. The source of systematic error differs significantly according to the 
method used to treat the heavy quark. In the Fermilab and NRQCD actions, the 
truncation of the perturbative and non-relativistic expansions is a major source 
of uncertainty. For the extrapolation method, the matching factor is known 
non-perturbatively and the non-relativistic expansion is not involved, but the 
extrapolation in the inverse heavy quark mass is a serious source of systematic 
error, as we discussed in Section 2. In addition, there is an error associated with 
the scale setting of the lattice, which is done by using some reference quantity — 
such as the rho meson mass, the pion (or kaon) decay constant, the string tension 
(or the "Sommer scale" ro which is also related to the heavy quark potential 
(93)), and the quarkonium mass spectra. In the quenched approximation the 
lattice scales determined from different quantities do not necessarily agree. This 
is a potential source of systematic uncertainty. 

More recently, there has been an attempt to control the ©[(amg)^] error by 
performing the simulations on a small physical volume, such as (0.4 fm)'^, and 
matching the results to larger lattices recursively (94). The result fs^ = 192 it 6 it 
4 MeV agrees with the previous calculations and the quoted error is significantly 
smaller than the present world average. Another attempt is a further refinement 
of the extrapolation method (95). By taking the continuum limit at each heavy 
quark mass, Rolf et al. eliminate the leading ©[(amg)^] error. The combined fit 
with the static result, which is also extrapolated to the continuum limit, yields a 
preliminary result /b^ = 206 it 10 MeV (95). By continuing in these and other 
directions, it will be possible to reduce the error in Jb to the 5% level. 

The B parameters have been calculated on quenched lattices via the static ac- 
tion (40, 96-98), the NRQCD action (42, 99, 100), and the extrapolation method 
(92, 101), for which the systematic error is better controlled by a combined fit 
with the static results (102). Their results are consistent with each other, and 
the current average is Bb = 1.33 ± 0.12 (12). 

The other B parameter Bs, which is relevant to the Bs meson width difference, 
has also been calculated in the quenched approximation (42, 100, 102-106). The 
physical result for AF^/F^ including next-to-leading order QCD corrections is 
0.074 lb 0.024 (107), which can be compared to the current experimental average 
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Omt'ofr (80). 

4-. 3 B^-B^ Mixing: Unquenching 

In the past few years, the unquenched calculation of fs has been performed by 
several groups (85, 86, 89, 108), and the results seemed larger than the quenched 
results by ~ 10%-15%. But the conclusion is now less clear because of an uncer- 
tainty in the chiral extrapolation as we discuss below. 

With the present algorithms to incorporate the effect of dynamical fermions the 
simulation cost increases as 1/m^, where ruq is the light quark mass (110). For 
example, for Wilson-type fermions the smallest possible sea quark mass is about 
a half of the strange quark mass rus even on the fastest available supercomputer. 
Hence, one must extrapolate the lattice data from the light quark mass region 
rUq > ms/2 toward the physical up and down quark masses, which are about 
m<j/25. In such an extrapolation the theoretical guide to restrict its functional 
form is given by chiral perturbation theory (ChPT) (111, 112), which provides a 
systematic framework to calculate physical quantities as an expansion around the 
massless limit with an expansion parameter oc nig. At next-to-leading (one- 
loop) order, ChPT predicts non-analytic behavior coming from pion loops. For 

instance, for the pion decay constant the non-analytic term is — ^^^jp-ln^ 
which is called the chiral logarithm. (/ is the pion decay constant in the massless 
limit.) Its coefficient depends only on the number of dynamically active quark 
flavors A'^^, which is a consequence of the chiral symmetry of QCD. 

The chiral logarithm could complicate the analysis of lattice data, because the 
effect of such non-analytic functional dependence can be missed if one simply 
relies on a polynomial fit of lattice data. Furthermore, when the lattice data are 
available only in the relatively heavy mass region niq > ms/2, the convergence 
of the chiral expansion itself is questionable. An explicit study using Nf = 2 
lattice data from Reference (113) suggests that the region niq > ms/2 is beyond 
the reach of ChPT (114). If this is true, one must introduce some model function 
to fit the lattice data, in order to make the chiral extrapolation consistent with 
ChPT, which introduces significant uncertainty in the chiral extrapolation. ^ 

The heavy mesons can be incorporated into the framework of ChPT. At next- 
to-leading order (121), the B meson decay constant depends on the light quark 
mass as follows: 

^f, _ 3(1 +3g^) ml .rnl , 

iw 4 — (i^TTF + + ^ 

/s 

Here, ^fg = fs V Mb and ^^j^ denotes its value in the massless light quark limit. 
The coupling g describes the B*B7r interaction in ChPT, which is experimentally 
measured through the D* Dtt decay as g = 0.59 ± 0.01 ± 0.07 (122, 123), if 
heavy quark symmetry is assumed. 

Figure 5 shows the chiral extrapolation of . The lattice data are from 
a recent unquenched simulation with two flavors of the 0(a)-improved Wilson 

^This point requires further investigation. Because the Wilson-type fermions exphcitly vi- 
olate the chiral symmetry, it is not clear a priori whether the lattice data can be fitted with 
the continuum ChPT formula. An extension of ChPT that includes such explicit chiral symme- 
try breaking has been proposed for the Wilson-type fermions (115-117) and for the staggered 
fermions (118-120). 
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Figure 5: Chiral extrapolation of ^/^^ {filled circles) and ^^/^^ (open squares). 
A polynomial fit is shown by solid lines; the fits with the hard cutoff chiral 
logarithm are shown for // = 300 {dotted curve), 500 {thin dashed curve) and oo 
{thick dashed curve) MeV. Quenched results are also shown {triangles). The plot 
is from Reference (124). 
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fermion (124). The NRQCD action is used for the heavy quark. In order to esti- 
mate the systematic uncertainty associated with the chiral logarithm, the lattice 
results are fitted with Equation 29, but the chiral logarithm term is modified as 
m^lnm^//x^ m^lnm^/(m^-|-^^). This functional form is motivated by a hard 
cutoff regularization of one-loop ChPT; it is designed such that the effect of the 
pion loop is suppressed above the cutoff scale /x, whereas Equation 29 is repro- 
duced when <C By taking a maximum parameter range /i = — > oo MeV, 
the results obtained are fs^ = 191(10)(+22) MeV, fs, = 215(9)(+}3) MeV, and 
Ibs/ fsa = 1-13(3) (1^2)) where the first error is statistical, and the second is the 
estimate of the uncertainty due to chiral extrapolation and other sources added 
in quadrature (124). Because the /^^ is hardly affected by the chiral logarithm, 
a large uncertainty remains only in /s^. The SU(3) breaking ratio fss/ fs^ con- 
tains significantly larger error than what was often reported in the quenched 
calculation, for which the effect of chiral logarithm was ignored. 

The problem of large uncertainty due to the chiral logarithm was also empha- 
sized in References (125, 126). A possible way to reduce the error is to extrapolate 
a ratio fnlf-K instead of the individual decay constants (127). Because the co- 
efficient of the chiral logarithm terms in and are numerically close, the 
uncertainty from the extrapolation will become smaller. For the cancellation of 
the chiral logarithms, a more natural choice is to consider fnlfu or the Grinstein 
ratio {fBs/fB)/{fDs/fD) (128), in which the chiral logarithm exactly cancels at 
the leading order of the 1/M expansion. These ratios are useful because the 
CLEO-c, an experimental progrom at Cornell, promises to measure and 
with a precision of a few per cent in a few years. An unquenched lattice calcula- 
tion of the Grinstein ratio has already been attempted aiming at a few per cent 
accuracy (129). 

To be convincing, however, the unquenched simulations must employ much 
smaller sea quark masses. With currently available computational resources, 
only staggered fermions are feasible as dynamical fermions to carry out such 
simulations. The first such calculation including the dynamical effects of one 
strange and two light quarks was presented recently for fs^ (130). The lightest 
sea quark mass is as small as ms/4 in this work. The result, = 260 it 29 MeV, 
is somewhat larger than the previous two- flavor result. 

The B parameter is also calculated in Reference (124). The chiral loga- 
rithm is less problematic here, because the coefficient of the one-loop chiral 
log term is -(1 - ^g^)/2 rather than 3(1 + ^g^)/'i in Equation 29 (131) and 
is numerically negligible in practice. Results for the physically relevant quan- 

titles are /b,/^ = 215(11)(1^^) MeV, /b./^ = 245(10)(1J?) MeV, and 
C = 1.14(3)(;i3). 

^.^ B D^*hu Form Factors 

The semi-leptonic decays b civ and h — > ulv may be used to determine the 
CKM matrix elements | Vcb \ and | Vub \ respectively, if their decay rates are theo- 
retically calculated in a model independent way. For the exclusive decay modes, 
such as -B — > D^*^lv and B nlv, the quantities to be calculated are their form 
factors. 

For the heavy-to-heavy decays, namely B D^*Hv^ one can use the heavy 
quark symmetry to restrict their form factors, i.e. several form factors with 
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different spin structures are related to one universal form factor called the Isgur- 
Wise function (132, 133). Furthermore, thanks to the symmetry between the 
initial and final states, the form factor is normalized to unity in the kinematical 
end point where the daughter D^*^ meson does not have spatial momentum in 
the rest frame of the initial B meson. Specifically, if one writes the differential 
decay rate in terms of w = v ■ v' {v and v' are the four-velocity of B and D* 
mesons, respectively) 

^nB_^ = m^>CMHu:f (30) 
aw 

with JC{w) a known kinematical factor, then the form factor J-{w) obeys J-{\) = 1. 
The CKM element \Vcb\ can thus be determined through the — > 1 limit of 
experimental data. 

However, because the heavy quark symmetry is exact only in the limit where 
both heavy quarks have infinite mass, one must take the 1/m corrections into 
account. Again using the heavy quark symmetry one can show that the leading 
correction to the relevant form factor /i^^(l), which is identical to Til) in the 
zero-recoil limit, is 0(l/m^) (134): 

/lAi(l) = ??a[1 + (^l/m.z + <Jl/m3], (31) 

with rjA a short distance correction that comes from the matching between QCD 
and HQET, and 

(2me)2 ^ 2m,2mb (2mb)2 ' ' 

The parameters ^p, iy, and Ia are non-perturbative quantities, which were pre- 
viously estimated using the non-relativistic quark model (135, 136). 

These three parameters appear in the 1/m expansions of other form factors 
because of the heavy quark symmetry (135, 137), and thus can be obtained from 
lattice calculation of double ratios of zero recoil matrix elements (21, 138) 

{D\c-f%\B){B\b-f^c\D) 



{D\c'y'^c\D){B\b-f^b\B) 
{D*\c'y%\B*){B*\h^c\D*) 
{D*\cj^c\D*){B*\h'^b\B*) 
{D*\c-f^j^b\B){B*\h^-f^c\D) 
{D*\c-f^j='c\D){B*m^j^b\B) 



7P[l-£pA + 0{l/m')]] , (33) 
7P[l-evA + 0{l/7n')]y , (34) 
{r?!f*[l-^AA + 0(lAn3)]}', (35) 



where A = (l/(2mc) — l/(2mfe)). The terms r/y* and 77^* are matching fac- 
tors between lattice theory and continuum HQET (139). The double ratios 
of matrix elements can be precisely calculated, because both statistical and 
systematic errors largely cancel in the ratio. In the quenched approximation, 
j^(l) = 0.913l|];[]^^l:[]:[j;^^ has been obtained (138). This may be compared to the 
quark model result (136) combined with a two-loop calculation of rj^ (140, 141), 
namely 0.907 ± 0.031, or the QCD sum rule result (142-144) 0.900 ± 0.038. 
The corresponding form factor G{1) for B Dlv is written as 



(36) 
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using the form factors /i+ and /i_ defined through 

{D{v')\cYb\B{v)) = VmBmD[{v + v')h+{w) + {v - v')h_{w)]. (37) 

The lattice calculation is more involved because it contains /i_(l), which cannot 
be obtained merely from the zero recoil matrix elements, but is possible using a 
similar technique (145). 

The Isgur-Wise function £,{w), which is a heavy quark limit of h^{w), can also 
be calculated on the lattice as a function of w. The lattice version of HQET that 
includes finite velocity was developed for this purpose (146) and was extended 
to include 1/m corrections (147-149). It is technically more difficult than the 
usual HQET, because the velocity receives finite renormalization, which must be 
determined from non-perturbative simulation (147, 150, 151). Lattice calcula- 
tions using the relativistic fermion have also been attempted (152, 153), but the 
control of systematic error at the level of a few percent (the precision necessary 
for the determination of \ Vcb\) is still challenging. 



4-5 B ^ Tilv form factors 



Because the symmetry between the initial and final states is lost in the heavy- 
to- light decays, the form factors are no longer normalized for B ttIv or plv. 
Therefore, lattice QCD must deal with the form factors themselves rather than 
just their 1/m expansion coefficients. 

The B —>■ ttIu form factors f^{q^) and f^{q^) are defined through 



{T,{k)\qrb\B{p))=f+{q^) 



{p + kf 



m. 



+ f{q' 



m 



B 



(38) 



where p and k are the momenta of the initial B and final vr mesons respectively, 
and = {p — k)'^ is a momentum transfer to the lepton pair. The minimum 
value of q'^, namely zero, occurs when the lepton and neutrino are parallel and 
the pion is energetically propelled to the opposite direction. The maximum value 



(ms 



occurs when the lepton and neutrino are emitted back-to- 



back with maximum energy and the pion stays at rest. The differential decay 
rate is written as 



dT{B -kIv) 
dq^ 



247r3 



\Vu,\\v.kf-mlflV\q')? 



(39) 



Because the discretization error becomes uncontrollable for momenta much 
larger than Aqcd ; the calculation of the form factors is feasible only in the large 
q^ (small recoil) region {q^ > 15 GeV^), where the spatial momentum of the pion 
is lower than roughly 1 GeV/c. The CKM matrix element \Vuh\ can be obtained 
by combining the experimental data integrated above some q^ value, say 15 GeV^, 
and the lattice results for the form factor |/'^(q^)P integrated in the same region 
with an appropriate kinematical factor. 

Extrapolation or interpolation in the heavy quark mass is done according to 
the heavy quark scaling law f~^{q^) ~ ^/rriB and f^{(f') ~ 1/ y/niB- The analysis 
can be more explicit if one works with the HQET motivated form factors fi{v ■ k) 
and f2{v ■ k). These are defined in the heavy quark limit as (154) 



{7rik)\qj^'b\B{v)) = 2 



fi{v-k)v>' + f2{v-k) 



ki" 
V ■ k 



(40) 
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Figure 6: q'^ dependence of B ^ -kIv form factors f'^{q^) {filled symbols) and 
f^{q^) {open symbols). Upward triangles, downward triangles, squares, circles 
denote data from the UKQCD (155), APE (156), Fermilab (157) and JLQCD 
(158) collaborations, respectively. 

where v'^ = p^^/ms is the four- velocity of the B meson. f^{q^) and f^{q^) are 
given by a linear combination of fi{v ■ k) and f2{v ■ k). 

Recently, five groups have carried out extensive quenched calculations of -B ^ 
ttIv form factors using the extrapolation method (155, 156), the Fermilab action 
(157), and the NRQCD action (158, 159) for heavy quarks. Results for the form 
factors f~^{q^) and f^{q^) are shown in Figure 6. Overall, results for f^{q^) 
are in agreement among four different groups within the error of order 20%, but 
some disagreement is observed for f^{q^). It originates from differences in the 
light quark mass dependence of the data. Further understanding is necessary to 
resolve these differences. 

The chiral extrapolation is more problematic for these form factors than for 
the decay constant, since the reference point q^ (or v ■ k) also varies with the 
light quark mass. If one sticks to a particular kinematical point, e.g., the zero 
recoil, the extrapolation must include a linear term in m-,^ in addition to the 
usual term. In fact, the soft pion relation f^{q^) = fs/ f-w is poorly satisfied 
unless the extrapolation includes the linear term (158). Furthermore, the 
chiral logarithm (160, 161) must be taken into account in future unquenched 
calculations. Some unquenched studies are in progress, especially using staggered 
dynamical quarks (162-164). 

4-6 B ^ plu and other form factors 

The decay B — > pliy can also be used to determine \Vub\- However, the lattice 
calculation of its form factors is more complicated than for B nliy, because 
there are four (not two) independent form factors corresponding to different spin- 
momentum combinations. Furthermore, the statistical signal of the Monte Carlo 
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Figure 7: Differential decay rate of the B — > pli' mode. Lattice data are shown 
with a fit curve discussed in the text. (From Reference (165) with permission). 



simulation is significantly worse for the rho meson, especially when finite spatial 
momentum is injected. A more fundamental problem is that the rho meson is 
an unstable particle that has a large width ~ 150 MeV. Lattice simulations of 
such an unstable particles should eventually treat the two-pion final states and 
take its phase shift into account, the practical feasibility of which is still an open 
question. 

Despite these problems, it is worth calculating those form factors on the lattice 
because the calculation provides a firmer theoretical guide for the experimental 
analysis than do the previous model calculations. A study was made by the 
UKQCD collaboration using the extrapolation method (165). (Earlier works 
on D —> K^*hu decays inculde References (166-169).) As in the B ttIu 
form factors the lattice calculation is reliable only in the high region. The 
extraction of \Vub\ can be made by combining a partially integrated decay rate 
with the lattice results in the corresponding region. The differential decay rate 
is written as 



dq^ 1927r3 



B 

(41) 

where A(g^) is a kinematical factor. H^{q^) and H^{q^) denote amplitudes with 
longitudinal and transverse rho meson polarization, respectively, and are written 
in terms of three form factors Ai{q^), A2{q^) and V{q^). Near the (/max o'^^ 
can parametrize these amplitudes as a^(l + h[q^ — ^max))- Figure 7 shows the 
lattice results for the differential decay rate with a fit curve using the above 
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parametrization. It demonstrates how lattice calculation can be used to compare 
theory with experiments. 

The systematic error is still as large as ~ 25% from the 1/M extrapolation 
and other sources. The calculation is in the quenched approximation and the 
chiral limit of the valence light quark is yet to be taken. These problems are 
being addressed. Indeed, a couple of large scale quenched lattice calculations are 
in progress using the extrapolation method (170, 171), and they have already 
reported preliminary results. 

Phenomenologically, the FCNC processes b — > s{d)^ and b s{d)l^l~ are also 
interesting because they occur through penguin and box diagrams, which are 
sensitive to possible new-physics contributions. Their exclusive decay modes are 
B K*j, B ^ pj, B ^ K^*H~^l~, and B K^*^^!?. For the radiative decays 
B — > K*^ and pj the lattice calculation of the form factor necessarily involves 
some model dependence, since the form factor at = is relevant. One must 
therefore extrapolate the lattice data from the high q'^ region assuming some 
functional form of the dependence, as was done in References (172, 173). 

For the B — > K^*h~^l~ and B K^*^uu decays, on the other hand, the lattice 
calculation may be useful in the large lepton invariant mass region. However, the 
long distance effect due to qq and cc resonances must be taken into account for 
the final states. 

4-7 HQET Parameters 

The HQET parameters, defined as 

f^liHQ) ^ -l—{HQ\QiiDfQ\HQ), (42) 
I^UHq) ^ -^{HQ\Qa-BQ\HQ), (43) 

appear in the analysis of the heavy quark expansion (174, 175).^ For instance, 
the inclusive decay rate of a heavy hadron Hq is written as 



1927r3 



( pI{Hq) - pI{Hq) \ , pUHq) ^ 
^ ' 2m^ / ^ m? 



Q 



(44) 

where the coefficients Cg and Cg are perturbatively calculable (176-179). Because 
the leading order term is independent of the type of hadrons, it gives interesting 
predictions of lifetime ratios: 



with a perturbative coefficient cg — 1-2 (180). 

The parameters p'^{Hq) and Pq{Hq) are non-perturbative parameters that 
depend on heavy hadrons Hq. Whereas the spin-dependent parameter Pq{Hq) 
is reliably determined from the hyperfine splitting ms* — rriB, the determination 
of p^{Hq) requires non-perturbative methods, such as lattice QCD. 

^Another notation, Ai = —jj.^ {B) and A2 — fJ-c (B) /3, is often used in the literature. 
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The lattice calculation of ^.'^{Hq) suffers from the subtraction of the quadratic 
divergence that appears in the matching of the lattice operator —QD^Q to its 
continuum counterpart. Because the perturbative subtraction involves large sys- 
tematic error (43), a non-perturbative subtraction has been attempted and the 
result Ai = 0.09 ± 0.14 GeV^ has been obtained (181, 182). Another possible 
approach is to rely on the mass formula, such as M§ = rrif, + A — Ai/2m;, for 
a spin-averaged IS B meson, and to fit the lattice data as a function of l/rub- 
The results of such analysis using the NRQCD action (183) and the Fermilab 
action (184) are Ai = -0.1 ± 0.4 GeV^ and -0.45 ± 0.12 GeV^, respectively 
In this method the quadratic divergence is subtracted away by the fitting, but 
the matching of the kinetic term in the non-relativistic expansion is a possible 
source of error. The reason for the disagreement with the above result with 
non-perturbative subtraction is not clear yet. 

It is also possible to avoid the subtraction of the quadratic divergence by con- 
centrating on the difference among different hadrons, e.g. iJ,'^{Af,) — /i^(-B), which 
is still useful in the estimate of lifetime ratios through Equation 45. A recent 
quenched calculation (185), /i^(Ab) — /i^(i?) = —0.01 ± 0.52 GeV^, is consistent 
with most phenomenological estimates. In particular, the well-known inconsis- 
tency of the heavy quark expansion of r(Ab)/r(i?) ~ 0.98 with the experimental 
value 0.77 ± 0.05 (80) continues to be a problem. At the l/m^ order, the spec- 
tator quark effect, which is expressed in terms of Ai? = four-quark operators, 
becomes important (180). A lattice calculation (186) of those matrix elements 
suggests that the spectator effects are indeed significant but not sufficiently large 
to account for the full discrepancy. 



4-8 B*Bn Coupling 

The heavy quark symmetry and chiral symmetry can be combined into the frame- 
work of heavy meson ChPT (187-189), which enables us to systematically expand 
physical amplitudes in terms of small pion momenta. This theory includes one 
additional parameter, g. It is related to the B*Btt coupling gB*Bn defined by 

{Bip)7T{q)\B*{p')) = -gB^Bnq^v''i27r)^S\p'-p-q), (46) 

where rj^ is the polarization vector of the B* meson, and the B*Btt coupling 
93* Bn is proportional to the low energy constant g as gB*B-K = '^gmriB/ f-w- The 
analogous coupling in charmed mesons can be measured experimentally from the 
D* Dtt decay width, yielding g = 0.59 ± 0.01 ± 0.07 (122, 123). Using this 
value in B physics is subject to uncertainties of order A/rUc. 

At tree level, the coupling is related to the semi-leptonic B — > vrZz^ decay form 
factor f~^{(f') in the large region. Using the soft pion theorem the form factor 
f~^{q^) near q^ ~ m^* behaves as 



Ib* gB*BTT 
2mB* (1 - q^/mB* 



Thus, a precise knowledge of gB*Bn can constrain the CKM matrix element \Vuh\- 
The coupling g also appears in the ChPT loop amplitudes as in the discussion of 
chiral extrapolation of Jb in Section 4.3. 

A direct lattice calculation of gB* b-k is possible if one modifies Equation 46 by 
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using the soft pion theorem as a matrix element of the axial-vector current 

{Bip)\q^A'^\B*{p',X)) = gs^Bn ^ J ^ml, (48) 

where q = p — p' ■ The lattice calculation of the matrix element on the left hand 
side is technically similar to the calculation of semi-leptonic form factors. A 
quenched calculation in the static limit was done several years ago and the result 
was g = 0.42 it 0.04 it 0.08 (190). More recently, Abada et al. performed an 
extensive calculation for the D*DTr coupling (191) and for B*Btt in the heavy 
quark mass limit (192). Their result g = 0.48 it 0.03 zt 0.11 (192) is consistent 
with the previous estimates. They also found that the l/mg dependence is not 
significant. 

4-9 Light- Cone Wave Function 

As discussed in Section 4.5, B decays with large recoil momentum are beyond 
the reach of current lattice calculations. Multi-hadron final states are even more 
difficult to treat on the lattice. 

For these energetic decays, a more natural theoretical treatment is to factor- 
ize the energetic interaction, which is perturbatively calculable, and the non- 
perturbative physics, which is governed by the energy scale of order Aqcd- Such 
a formalism was developed by Beneke et al. a few years ago (193, 194). For 
example, the amplitude of the non-leptonic B vrvr decay is written as 

{7r{p')7riq)mBip)) = f^'^^q^) C duTl {u)^^{u) 

Jo 

+ / dCdudvT,"{C,u,v)<^B{0'^Au)'^Av)- (49) 
Jo 

Here Qi is a four-fermion operator and T's are the perturbatively calculable short- 
distance part, whereas f^^'^ is the B — > ttIu semi-leptonic decay form factor, 
and $'s are the light-cone wave functions. Thus, the genuinely non-perturbative 
quantities are these light-cone wave functions. 

There are lattice studies on the light-cone wave function of the pion. The light- 
cone wave function (f)-,^{u,Q") is defined by the inverse Fourier transform of the 
matrix element 

(0|J(0)e/°'^^-^W7A.75^x(x)|7r+(p))|p2=o = -iP^.U Cdue-'P^(t)^{u,Q^). (50) 

Jo 

The n-th moment of 4>tt{u, Q"^), defined as {^^) = d^^'^4)T^{^, Q^), can be related 
to the pionic matrix element with the local bilinear operator with n derivatives 
as follows: 

(0|d7/.75 D^,,■■■ Df,„ u\tt+{p)) = i''Up^,p^,, • • (51) 

The second moment of the light-cone wave function has been calculated by several 
lattice groups (195-200). A similar calculation is needed for the B meson in order 
to establish contact with the B physics phenomenology. 

A more direct calculation of the light-cone wave function, not limited to its 
moments, has also been proposed (201).^ A pioneering numerical study has been 
carried out using this method (203), which is worthy of further studies. 

■^This method can also be applied to the calculation of the B meson shape function (202), 
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5 Future perspectives 

Lattice QCD has many applications, but the flavor physics is where it plays a 
crucial role in the exploration of physics beyond the Standard Model through 
model independent calculation of weak matrix elements. Because the precision 
of the calculation is crucial in such studies, the systematic errors in the lattice 
calculation have to be well understood and reduced as much as possible. In the 
past decade many important steps have been made in this direction: 

• O (a) -improvement of lattice actions and operators. The discretization er- 
rors are described by Symanzik's effective theory, which also provides a 
method to improve lattice actions in a systematic way (204). It is now 
common to use the 0(a)-improved action, for which the leading discretiza- 
tion error is O(a^). Further improvement is also being pursued. 

• Effective theory for heavy quarks. The ideas of HQET and NRQCD are 
naturally implemented on the lattice. The problem of dealing with heavy 
quarks on the lattice has essentially been solved (15). 

• Improved perturbation theory. Perturbation theory is still necessary in the 
lattice calculation to match the lattice action onto the target continuum 
theory. The convergence of the perturbative expansion was dramatically 
improved by the use of the renormalized coupling (205). 

The first two advances are related to the separation of different energy scales by 
the use of effective theories, and the third is needed for the matching of these 
effective theories. 

We have reviewed the results that have relied on these methodological devel- 
opments. One of the main challenges now is to further substantially reduce the 
systematic error. This could be achieved by pushing the idea of the effective 
theories, i.e. including the higher order improvement terms and higher order 
relativistic corrections. These have to be done together with the matching with 
the higher order perturbation theory. Although the higher order perturbative 
calculation is technically rather demanding on the lattice, some work is already 
in progress (206). Another proposed way to reduce the systematic errors, which 
does not rely on perturbative matching, is the recursive non-perturbative match- 
ing methods (44, 94). 

To obtain truly model independent results from lattice QCD, the quenched 
approximation must be abandoned. Large scale simulations that include the 
effects of dynamical quarks have already been carried out by several lattice groups 
and will be performed more extensively in the future. One of the important 
effects of dynamical quarks is the pion loop effect, which leads to the chiral 
logarithm and affects the chiral extrapolation of lattice data in a non-trivial 
way. As a result, the systematic uncertainty due to the chiral extrapolation is 
still quite significant. Although our discussion on this problem used fs as an 
example, the same is true for almost all other lattice observables as long as their 
chiral behavior is known from ChPT. (Otherwise, there is no theoretical guide 
for the functional form of chiral extrapolation.) To eliminate the uncertainty 
from the chiral extrapolation, one must push the dynamical quark masses in 
the lattice simulations down to the region where the chiral logarithm becomes 

which describes the effect of the Fermi motion of b quark on the inclusinve B decays B Xsf 
and B — > X^lv- 
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a visible effect. It is presumably lower than ms/3. Such simulation is currently 
feasible only with the staggered fermion action for dynamical quarks, which is 
numerically so cheap that the MILC collaboration achieved even m^/G. The 
results for several physical quantities agree with their experimental values quite 
remarkably (58). The price, however, is the introduction of fictitious species, 
sometimes called "tastes." The taste breaking effect introduces an additional 
source of systematic uncertainty. Furthermore, in the dynamical simulation the 
fourth root of the fermionic determinant must be taken, and it is not known 
whether this gives a local field theory in the continuum limit (207). These points 
are still controversial. Other lattice fermion formulations should also be pursued, 
for which the developments of simulation algorithms are essential. 

As we have discussed, in order to achieve the goal of precise, model-independent 
simulation of heavy quarks, improvements are necessary throughout lattice QCD, 
i.e. from the higher order perturbative calculation to the simulation of really light 
dynamical quarks. Despite the difficulties, these future improvements will be 
worth the effort, because from an improved flavor physics one may probe physics 
beyond the Standard Model. 

The application of lattice QCD to heavy quark physics is not limited to /b 
and Bb- Some other relevant quantities, such as semi-leptonic decay form fac- 
tors and heavy quark expansion parameters, have already been studied on the 
lattice. There are, however, many other important quantities that require non- 
perturbative calculations. The FCNC processes are sensitive to 

new physics, and hence the calculation of their form factors has direct relevance 
to the new-physics search. The two-body decays of the B meson have become the- 
oretically tractable since the factorization of short and long distance physics was 
proven. Non-perturbative inputs are necessary for the long distance part, which 
is the light-cone distribution function. Inclusive B decays, such as B ^ Xg'j and 
B — > Xuliy, can be calculated perturbatively by using the heavy quark expansion, 
but in the kinematical region of interest the simple OPE breaks down, so that 
resummation of a certain class of higher-twist terms is necessary. Thus, a non- 
perturbative calculation of the quantity called the shape function, which describes 
the Fermi motion of the bottom quark, is required to predict their energy dis- 
tributions. Aside from the B decays, the heavy meson and baryon spectroscopy 
is also attracting attention, since the recent discovery of such exotic particles as 
Dsj{2317) and X(3872). Lattice QCD is the prime tool for the non-perturbative 
study of these particles. All these applications will be investigated in the future. 
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